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We present first-principle numerical calculations for few particle solutions of the attractive Bose- 
Hubbard model with periodic boundary conditions. We show that the low-energy many-body states 
found by numerical diagonalization can be written as translational superposition states of compact 
composite systems of particles. These compact states break the translational symmetry of the 
problem and their center-of-mass and internal excitations offer simple explanations of the energy 
spectrum of the full model. 



I. INTRODUCTION 

Interacting many-body quantum systems present theo- 
retical challenges in all branches of physics and chemistry, 
and a wide range of methods have been developed to pro- 
vide efficient computational methods and to gain insight 
in their detailed structure and dynamics. Separation of 
different physical degrees of freedom is a widely used ap- 
proach, e.g., in chemistry, where the Born-Oppenheimer 
approximation treats the electronic motion in the pres- 
ence of classically fixed nuclei, and subsequently the re- 
sulting energy levels together with the direct nuclear 
interaction constitute potential energy surfaces for the 
quantum motion of the nuclei [Ij. The nuclear motion, 
in turn, can be separated into vibrational motion around 
minimum energy configurations, and quantum transla- 
tional and rotational motion of the entire molecule. The 
latter reflect the translational and rotational invariance 
of the Hamiltonian and establish total momentum and 
angular momentum as conserved quantities with associ- 
ated good quantum numbers. 

A similar separation of the interaction energy asso- 
ciated with the relative motion and with the transla- 
tion and rotation of the composite object of particles is 
used for the state of atomic nuclei [2j, and it has also 
been applied in the interpretation of excited state spectra 
of atoms with spatially (anti-) correlated electrons [SHS]. 
The validity of such a separation in relative motion of 
individual particles and collective motion of the bound 
composite object is usually explained with reference to a 
separation in energy scales which suppresses the coupling 
between the different degrees of freedom and justifies an 
adiabatic separation, e.g., of the slow nuclear and fast 
electronic motion in molecules. 

In ultra-cold atom physics it is possible to study the 
case of binding due to very weak and tunable interac- 
tions, and we may probe the transition between a collec- 
tion of independent particles and a single massive object 
held together by the attractive interactions. Attractive 
bosons may thus form a so-called bright-soliton, which in 
mean field theory [6]-[8] is described by a localized solu- 
tion to the non-linear Schrodinger Equation. The choice 
of a particular localization breaks the translational sym- 
metry of the system in free space, but rather than restor- 
ing translational symmetry by diluting the atoms in a 



uniform mean field wave function, a quantum superposi- 
tion state of translationally displaced copies of the soliton 
wave function may form a better Ansatz for the many- 
body quantum state. 

A similar situation occurring for atoms in a trap was 
analyzed by Pethick and Pitaevskii [9j, see also p!QHT2] . 
where again a compact bright soliton like object may 
form, and where this object may exercise quantum mo- 
tion described by the center-of-mass ground state wave 
function in the trap. If the quantum mechanical excur- 
sions of the center-of-mass motion are comparable to or 
larger than the soliton width, the state is poorly de- 
scribed by a single mean field wave function, and as a 
consequence the high degree of order in the system is not 
reflected by a macroscopic population of a single particle 
state. 




Figure 1: The Bose-Hubbard model describing atoms popu- 
lating single particle states localized at lattice potential min- 
ima. Upper panel: The atoms experience tunneling between 
sites, while atoms on the same site interact with interaction 
strength [/ < 0. Lower panel: Many-body quantum state in 
a finite lattice potential with periodic boundary conditions. 
The translated copies of a classical configuration of atoms 
depicted as red balls illustrate how the full state is a super- 
position of translated many-body quantum states. Similarly, 
the green curve shows the atomic density in a localized state, 
while the density in the translated replica are indicated by 
the blue curves. 
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In [131 m mean-field and Bogoliubov theories are used 
to reveal how a system of bosons on a ring with periodic 
boundary conditions breaks the translational symmetry 
and undergoes a quantum phase transition towards a lo- 
calized soliton-like state, and how critical fluctuations 
and localization may occur due to a weak symmetry- 
breaking perturbation or measurements on the system. 
Anderson localization of the mean field solution in an ex- 
ternal disordered potential on length scales much larger 
than the soliton width is reported in \T5]. 

The separation between relative degrees of freedom and 
center-of-mass motion of composite, bound objects, does 
not rely on mean field theory. In [16j it is thus shown that 
the soliton-like object can scatter and form macroscopic 
superposition states of transmitted and reflected compo- 
nents, while a many body description is used to describe 
the decay of the soliton into fragmented condensates in 
[T7]. Two particle bound states on a lattice separate ex- 
actly in center-of-mass and relative coordinates [18-20J, 
and in [21], it is shown how the ground state solution to 
the Lieb-Liniger model of an attractive Bose gas in one 
dimension can be precisely represented as a translational 
invariant state with a definite dependence on the rela- 
tive coordinates between the particles following a Bethe 
Ansatz with exponentially decreasing terms. 

In this article, we will analyze the attractive Bose- 
Hubbard model on a finite lattice with periodic bound- 
ary conditions as illustrated in Figure [l] We will treat 
the case of few atoms and few sites, which permits ex- 
act numerical diagonalization of the problem, and we will 
subsequently analyze the results, and show how the dis- 
crete translational symmetry of the problem gives rise 
to states, which can be represented as superpositions of 
translated versions of compact symmetry breaking states. 
We will study both the ground state and excited states 
of the system, and we will identify how the low-energy 
part of the energy spectrum can be understood in terms 
of the center-of-mass and relative motion of the atoms. 

In Sec. [TT| we introduce the system and identify the 
symmetries of the problem which enable a reduction of 
the full many-body problem to a numerically tractable 
form. In Sec. |III[ we present numerically determined en- 
ergy spectra and correlation functions in support of our 
separation of the problem. In Sec. |IV[ we show that 
our states can be represented by the desired Ansatz. In 
Sees. |V| and \VJ\ we show how the underlying compact 
quantum states are related to and can explain features 
in the computed energy spectrum for the model. We 
conclude in Sec. IVIII 



II. THE SYSTEM AND ITS SYMMETRIES 

We consider a one-dimensional lattice with M sites 
numbered 0,1,2,...,M — 1, and we restrict our anal- 
ysis to the case where atoms can access only one single- 
particle state at each site j. A system of identical 
bosons on such a lattice is most conveniently treated in 



second quantization, and in the tight-binding approxima- 
tion, the system is well-described by the Bose-Hubbard 
Hamiltonian 

H = -JY^ + H.c.) + ^^fijihj - 1), (1) 



where J and U are the tunneling and the on-site inter- 
action strengths, respectively. In the present paper we 
consider attractive interactions, U < 0. The operator 
Oj creates a particle in the single particle state and 

fij = a^jCLj is the corresponding number operator on the 
jth site. Periodic boundary conditions are ensured by 
letting all site indices be implicitly understood as mod- 
ulo M. 

H conserves the number of particles, and we will per- 
form our calculations within fixed Hilbert spaces, , 
of states with TV particles on M sites, conveniently 
spanned by the basis of multi-mode Fock states 

• + riM-i = N. There are 



|no, . . . ,nM-i) where no 



■ N 



N 



(2) 



elements in the set of such basis vectors, and in Ap- 
pendix [Aj we discuss an ordering of these basis states 
which makes the action of the Hamiltonian particularly 
easy to calculate. 

The unitary many-body translation operator, T, acts 
on the many-body state by moving all particles one site 
"to the right": 

T |no, . . . , riM-i) I^M-i, ^0, ^i, • • • , , (3) 

and when we apply it M times on any of the Fock basis 
states, we return to the original state, i.e. 



= 1. 



(4) 



The action of the translation operator can also be spec- 
ified by the relations 



Ta] = aj^iT, 



(5) 



and their hermit ian conjugates. It follows directly that 
the Bose-Hubbard Hamiltonian commutes with T and 



thus simultaneous ei 



genstates exist for the two operators. 



Since the relation (|4) dictates that the T-eigenvalues are 
Mth roots of unity, we can label those eigenstates 
as 



(6) 



In the limit of a large number of sites, the energies 
form bands labeled by the excitation number m = 1, 2, . . . 
and by the quasi-momentum 



Qn 



27rn 



g] - 7r;7r], 



(7) 
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(a) UjJ = (b) U/J = -0.25 (c) U/J = -0.5 (d) U/J = -1 




Qn Qn Qn Qn 



Figure 2: (Color online) Energy spectra for 11 particles on an 13 site lattice are shown as functions of the quasi-momentum qn 
for (a) [//J = 0, (b) U/J= —0.25, (c) U/J= —0.5, and (d) U/J= —1. The colored solid lines connect the energies computed 
at the discrete allowed values of qn- 



which, acquires M discrete values in a finite lattice. 

Finally, the system is invariant under parity symmetry, 
which ensures that = and which imposes con- 
straints on the real and imaginary parts of our expansion 
coefficients that are useful in the following. 

We can divide the basis Bf^ into P disjoint equivalence 
classes of basis vectors 



we can choose such a basis as 



i5j^ = f«u£:(2)u---uf(^) 



(8) 



where and lie in the same equivalence class 

£^^^ iff = for some k. Every state |^) can 

be expanded uniquely as 



p \s^^^\-i 



(9) 



where we have renamed the basis states such that 

£:(^-) = {|^^^-)),|^l^-)),...,|$g^^^ (10) 



\£^^^\ being the number of elements in S^^\ and \^ 

If 1^) is a translation eigenstate with quasi-momentum 
Qn, then for any two elements |^a^) and l*^*^"^^) from 

and the expansion (|9| takes the simpler form 

p \s^'^\-i 



j=l fe=0 



(11) 



where Cj accounts for the overlap between |^) and states 
from S'^^^ 

Since H and T commute, H is block diagonal in a basis 
of translation eigenstates. Inspired by the expansion ( 11 ) 



M 



E ^ 

k=0 



-27ri/cn/M 



(12) 



with n = 0, . . . , M — 1 and j = 1, . . . , P. These states are 
orthonormal in both indices and the matrix representa- 
tion of H is block-diagonal. The problem thus reduces 
to the identification of eigenvalues and vectors of M ma- 
trices of dimension Vf^ /M rather than of one matrix of 
dimension Tj^. 



III. NUMERICAL RESULTS 

Using the state vector representation described in the 
Appendix, and the simplification of the problem following 
from the above analysis, the system becomes amenable 
to standard numerical diasonalization routines, and the 
eigenstates and -energies ([6| of ^ can be found. 



A. Energy spectra 

In Figure |2| the energy spectra for 11 particles 
on a 13 site lattice are shown for different interac- 
tion strengths delineating the transition from the non- 
interacting regime to the strongly attractive regime. The 
colored solid lines are used to connect the energies at the 
discrete allowed values of Qn in the figure. 

In part (a) of the figure, U/ J = and the interac- 
tion term of H vanishes. Thus, the system can be solved 
by finding the single-particle eigenstates and -energies 
on the discrete lattice and forming A/'-particle product 
states. Weak interactions among the atoms are intro- 
duced in part (b) of the figure. They couple the product 
states, and as the interactions get stronger in parts (c) 
and (d) of the figure, the spectrum develops structures 
at different scales. These are the structures that we will 
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associate with different modes of atomic motion: low en- 
ergy translational motion of the entire composite object, 
shghtly higher energy motion of atoms outside the bound 
composite object, and large energy gaps associated with 
the excitation of different numbers of atoms out of the 
bound composite object. 

In the subsequent sections, we examine this separation 
hypothesis further by an analysis of the eigenstates found 
by our diagonalization, and we shall also quantitatively 
address the energy contributions from center-of-mass and 
individual particle motion. 

Here, let us consider the coarse scale structures in Fig- 
ure [2jd) The interaction energy of n particles on a single 
site is ^n(n — 1), and this explains the coarse separation 
of the spectrum into bands of states around ^N{N — 1), 
f (AT - 1 ) (TV - 2) , and f (TV - 2) (TV - 3) , which are transla- 
tional superpositions of states where either all atoms oc- 
cupy the localized condensate state, or one or two atoms 
are excited out of the condensate. In the total absence 
of tunneling, U/J — oo, each of these bands become 
degenerate, while in Figure [2](d), in the presence of both 
tunneling and strong attraction, the degeneracy is lifted 
within each quasi- momentum subspace, and the bands 
split into a series of almost flat sub-bands. The crossing 
structures observed at the weaker interaction strength in 
Figure |2|c), will be discussed below. 



B. Eigenstates 

All the energy eigenstates found by our diagonaliza- 
tion are quasi-momentum eigenstates, and hence the 
particle density is uniformly distributed over all lat- 
tice sites. Usually, the structure within the many- 
body states is analyzed by correlation functions, and 
in Figure [sj we show the flrst order (amplitude) coher- 
ence function, (^5^|a^a/c|^5^), (the one-body density ma- 
trix), and the second order (density) coherence function, 
(^J^Iaja^a/cttjl^J^). The functions are shown for the two 
lowest energy states with vanishing quasi-momentum. 
The upper plots reveal the constant density along the 
identical j = k diagonal elements, and they show the van- 
ishing of the flrst order coherence function for large sep- 
aration. Diagonalization of the one-body density matrix, 
yields M almost identical eigenvalues around N/M^ indi- 
cating that the atoms are evenly distributed on M single 
particle modes, and according to the Penrose- Onsager 
criterion, the system of bosons does not populate a con- 
densate. The width of the diagonal bands in the flgure 
indicate, however, that these single particle modes are 
localized on a small number of neighboring sites. This is 
similar to the observation by Pethick and Pitaevskii [9j 
that trapped attractive gases do not condense according 
to the conventional criteria of Bose-Einstein condensa- 
tion. 

The lower plots show the density coherence functions 
with the intuitive interpretation as the probability of 



6 



2 

^ 

-2 
-4 
-6 
6 



2 

^ 

-2 
-4 
-6 



j j 




0.5 1 1.5 2 



Figure 3: (Color online) First and second order coher- 
ence functions, (^^^^ |a]afc ) and {'^'^\a^jalakaj\'^'^) , for 
the two lowest many-body eigenstates with vanishing quasi- 
momentum for U/ J — —0.25. 



flnding a particle at a given site, /c, given that another 
particle has been seen at another site, j. The large values 
along the diagonal show that, while evenly distributed 
on average, the atoms are not independent, but cluster 
strongly around each other within a range comparable to 
the single-particle coherence range in the upper panels. 
This of course supports our hypothesis, that the atoms 
form a compact condensate wave function, like a bright 
solitonic wave in mean fleld theory, but in order to min- 
imize its kinetic energy, the center-of-mass of this wave 
function should be completely delocalized over the lat- 
tice. 

We note, that localized bright or dark solitons in mean 
fleld theory break the translational symmetry of the un- 
derlying equations, while Bogoliubov theory returns both 
quasi-particle excitation modes and two Goldstone modes 
associated with the U{1) phase symmetry breaking and 
the breaking of translational symmetry. The latter Gold- 
stone mode causes a growing position uncertainty of the 
whole mean fleld solution and shows explicitly that local- 
ized solitons are not stationary solutions [22, 23j. Simi- 
larly, in our case, a localized many-body state would not 
be an energy eigenstate. We do not have recourse to Bo- 
goliubov theory, but in the following section, we shall see 
how the same conceptual idea works in our exact solu- 
tions of the Bose-Hubbard model. 
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Figure 4: (Color online) Complex center-of-mass coordinates 
of all the Fock basis states for 11 modes and 4 particles. The 
gray wedge encapsulates states with a center of mass coordi- 
nate with a complex phase in the interval ]— tt/M; tt/M]. 



IV. DECOUPLING OF THE TRANSLATIONAL 
MOTION OF EIGENSTATES 

The goal of this section is to identify localized, compact 
states of the N atoms, which enable us to write the 
translation eigenstates \^^) as superposition states of 
the form. 



M-l 



(13) 



/c=0 



where C is a normalization constant. The lower part of 
Figure [T] attempts to illustrate the idea behind an Ansatz 
of precisely this form, where the phase factors in ( 13 ) lead 



to a translation eigenstate with quasi- momentum q^. 



A. Superposition and projection 



We introduce Eq. (13) as a superposition of translated 
replica of the state but we may also read the equa- 

tion in a different manner by noticing that the self-adjoint 
operator 



M-l 

M ^ ' 

k=0 



(14) 



is a projection operator, which projects any state in Hj^ 
onto the subspace of translation eigenstates with quasi- 
momentum g^. The superposition (13) of translated 



replica of a single localized state is thus, up to normal- 
ization, identical to the projection of that state on its Qn 
quasi-momentum eigenstate component. 



An immediate consequence of this is that given the tar- 
get state 1^5^) there are infinitely many possible choices 
of |q^5t^), and all the states forming, e.g., the lowest energy 
band in Figure |2] can be constructed from the same state 
in a quite tautological manner, by first adding all these 
states together in a superposition with arbitrary non- 
vanishing expansion coefficients, and subsequently apply- 
ing the operators (14) to extract all the eigenstates one 



by one. This is merely a formal property of states obey- 
ing a symmetry, and hence it does not provide a physical 
explanation of the spatial correlations in the states. 

Rather than forming the underlying state from 
the energy eigenstates in each excitation band, we will 
build their formation on the physical principle, that they 
are well-localized on the lattice and centered at the site 
j = 0. Since each of the equivalence classes introduced 
in (|8| contains states which are related to each other by 
displacement, they all project onto the same state under 
the operation by f^, which in turn implies, that we can 
build any quasi-momentum eigenstate from (13) 

with composed of only one member from each equiv- 
alence class 



lC)>€f(^). 



We now select the representatives of each equivalence 
class according to their spatial location on the lattice. 
Due to the periodicity of the system we can picture the 
sites lying in a circle as in Figure [T] and thus assign a 
two-dimensional position (xj^yj) to each of the M sites. 
For convenience we use complex notation Zj = Xj + iyj = 
g27rij/M j^Q^ define a center-of-mass coordinate 

as the arithmetic mean of the coordinates in the com- 
plex plane. This definition of the center of mass has the 
advantage, that when we apply the translation opera- 
tor on a basis state, we rotate the complex phase an- 
gle of all particles by 27r/M radians counter-clockwise — 
and the center of mass thus rotates by the same angle 
in the complex plane. This rotation, couples the dif- 
ferent elements within each equivalence class (10), and 



we propose to expand our candidate localized state on 
members from each equivalence class which have their 
center-of-mass complex phase angle within the interval 
] — tt/M; tt/M]. Figure [4] depicts the center-of-mass co- 
ordinate for all Fock states with 11 sites and 4 particles, 
and the shaded wedge indicates the values for which we 
include the state in our expansion. In practice for the 
state vectors found by our numerical diagonaliza- 

tion, we merely retain the expansion coefficients on the 
basis states with center-of-mass in the wedge. To sum- 
marize, we choose the underlying state as 



l"n) 



(15) 



where W is the projection onto the space spanned by the 
Fock states with center-of-mass within the wedge, and 
even superposition states of those Fock state pairs with 



C//J = 



U/J = -0.1 



U/J -0.25 
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Figure 5: (Color online) Densities of \a^) for the two lowest energy excitations {m =1,2) for four different attractive interaction 
strengths. The rightmost panels correspond to the two lowest excitations in Figure [2|c). 



cent er-of- mass on the upper and lower angular limits of 
the wedge. 



V. LOCALIZED STATES 

In the previous section we motivated a procedure 
to find localized states, which yields the energy 

eigenstate, by the translation and superposition 



Ansatz (13) 
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Figure 6: (Color online) Energy levels and behavior of local- 
ized states near avoided crossings. The upper panel shows the 
energy levels of the 2nd to 5th excitations dXU / J — —0.45. In 
the middle and lower panels we show the overlap of the local- 
ized states identified for each quasi-momentum with the one 
for the definite value no = —3 in the 4th and 3rd excitations, 
respectively, marked with circles in the upper panel. 



lO 



(16) 



Figure [5] shows the densities of these for four different 
interaction strengths and for the two lowest energy exci- 
tation numbers, m = 1,2. For low interaction strength, 
tunneling is significant, which causes the states to spread 
out, while for larger attraction the states narrow to min- 
imize energy, and in the limit U/J — oo all particles 
are located on one site in the m = 1 states. The figure 
also shows, that even though the states |q^5^) are deter- 
mined independently for each state I^J^), they are very 
similar for different g^, supporting our assumption that 
the center-of-mass and the internal relative motion of the 
atoms separate. 

In our previous discussion we justified the separation 
of different motional degrees of freedom by their different 
energy scales. Parts (a) and (b) of Figure |2]lend no sup- 
port for the separation, while parts (c) and (d) display 
a clear distinction between the energy scales associated 
with the center-of-mass motion and the relative motion. 
A magnified view of excitations 2-5 in Figure [2|c) are 
reproduced in the upper panel of Figure |6] (for a slightly 
different interaction strength). Among these states an 
avoided crossing is observed between the 3rd and 4th ex- 
citation bands at zero quasi-momentum. This suggests 
that the separation Ansatz does not hold here, and in- 
deed, the localized states identified for these states vary, 
as we pass the crossing region. 

This is quantified in the two lower panels of Figure [6j 
where the overlap of the localized states I^J^) obtained 
from different quasi-momentum states with a fixed local- 
ized state \a^^) is plotted. The solid blue curve in the 
middle panel shows that the localized states are almost 
constant for negative quasi-momentum states in the 4th 
excitation band with a near unit overlap with the typical 
state lo^is) in that interval. However, the states abruptly 
changes character at the crossing, and acquires the char- 
acter of the negative quasi-momentum states in the 3rd 
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excitation band as seen from the dashed blue curve in 
the lower panel. Consistently, the same behavior is ob- 
served for the states in the 3rd excitation band. This 
is a typical example of the diabatic transition between 
pairs of states near degeneracies in a physical system, 
and we may in fact imagine an experimental application 
of this crossing phenomenon, where the system is driven 
through the avoided crossing by acceleration of the lat- 
tice, and thus superpositions of two different many-body 
quantum states may be controllably prepared. 

The curves in the upper panel suggests that further 
avoided crossing with the 2nd and 5th excitation bands 
occur at the band edges — but more pronounced in the 
spectra for weaker attraction — and the overlap ampli- 
tudes in the lower panels confirm a non-adiabatic transi- 
tion behavior toward these states. 



A. Many-body character of the localized states 

We have shown, that the energy eigenstates in the 
Bose-Hubbard model in the case of strong attractive in- 
teractions can be written as translational superposition 
states of localized composite many-body states of the 
atoms. Let us now investigate their many-body char- 
acters in more detail. 

For each state we can compute the one-body den- 
sity matrix p with the matrix elements 



(17) 



This quantity is equivalent to the first order coherence 
function, but unlike the results shown in Figure [3) here 
we compute it for the localized state. In Figure [7| is plot- 
ted the three largest eigenvalues of p for the states \o^^) 
with n = and m corresponding to the ground band 
and to the lowest sub-band of the first and second ex- 
cited band. We we see that the largest eigenvalue of p 
is comparable to the number of particles, which implies 
that in each state many of the atoms occupy the 

same single particle wave function, which in turn is given 
by the corresponding eigenvector of p. In the limit of 
strong attraction U/J — oo, for the ground state the 
largest eigenvalue Ai approaches the number of particles 
= 11, while for the states in the first and second ex- 
cited bands, the largest eigenvalue approaches TV — 1 and 
A/" — 2, respectively. This is a clear indication, that in the 
ground band, \a\) is an A^-particle condensate with all 
particles in some single particle state \Xn)i whose density 
is easily extracted from Figure [5j while in the first and 
second excited band, one and two atoms are excited out 
of the condensate, respectively. This is similar to the un- 
derstanding offered by the (number-conserving) Bogoli- 
ubov treatment of interacting bosons [24J in the limit of 
particle-like excitations out of the condensate wave func- 
tion. 

The middle panel of Figure [7| shows that in the first ex- 
cited band, the second largest eigenvalue of p is unity for 
a wide range of strong attraction, suggesting that a single 



excited atom occupies a definite single particle state \cr^) 
orthogonal to the condensate mode |x5?)- Qualitatively, 
this particle tunnels among M lattice sites, exposed to 
an attractive potential dip at the site occupied by the 
condensate, and we expect the particle to have discrete 
eigenstates equivalent to plane waves with a suppressed 
density at the condensate location. The condensate wave 
function and this excited state wave function are indeed 
available as the eigenvectors of p corresponding to the 
eigenvalues Ai and A2, respectively, and in Figure [S) we 
show these functions for U/ J = —1.75 for the three lowest 
excitations in the first excited energy band (m = 2, 3, 4). 
The condensate wave function shown in the upper panel 
is strongly peaked with small wings due to the tunnel- 
ing, and the singly populated wave functions for the three 
lowest states of the first excited band are indeed of the 
expected form with an increasing number of nodes corre- 
sponding to scattering states on the ring in the presence 
of the condensate. 

For states in the second band the largest eigenvalue of 
p has the limiting value N — 2^ but the two subsequent 
eigenvalues show a more complicated behavior. This is 
compatible with a condensate state with N — 2 atoms, 
while the two excited atoms populate a non-trivial two- 
particle state. To investigate this state further, be is ap- 
plied N — 2 times to |q^5^), where be is the annihilation 
operator that removes a particle in the condensate mode 
Ix^). The wave function of the resulting two-particle en- 
tangled state is illustrated in Figure [9] for the three 
lowest zero quasi-momentum states in the second excited 
manifold. The states are shown for two different inter ac- 





Second excited band 
First excited band 
Ground band 




U/J 



Figure 7: (Color online) The three largest eigenvalues (Ai > 
A2 > A3) of the one-body density matrix (17) for the localized 



states \a^) with the lowest energy in the ground band, the 
first and the second energy band. The limiting value of Ai for 
strong attraction equals the number of particles, N = 11, in 
for the ground band. The rich structures in the eigenvalues at 
weak interactions are due multi-level intraband level crossings 
which change the character of the lowest eigenstates in the 
first and second energy bands. 
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U/J = -1 



U/J = -2.5 



Figure 8: Wave functions of the single-particle states found as 
eigenvectors of the one-body density matrix (17) for \a^) with 
n = and different values of m and with U — —1.75. Similar 
results are found for all values of n, but not with real- valued 
wave functions. The upper panel shows the condensate wave- 
function which are practically indistinguishable for different 
values of m in the first band, and in the lower panel we see 
the wave function of the single excited particle in the three 
lowest states in the first energy band. 



tion strengths and are evidently not (symmetrized) prod- 
uct states and can therefore not be regarded as two par- 
ticles occupying two (different) single-particle states. In- 
deed, the attractive interaction strongly correlates the 
atoms on the same side of the condensate peak at j = 0, 
and as the attraction increases the wave-function be- 
comes more and more narrow around the diagonal in the 
plot, indicating an increased spatial bunching of the two 
particles. 

We note that the conventional Bogoliubov treatment is 
formally an expansion in 1/ y/N [24j , and interactions be- 
tween atoms outside the condensate wavefunctions can be 
neglected in the limit of large N . Due to the near degen- 
eracy of non-condensate single particle states shown in 
Figure |8) these interactions are significant in our problem 
and cause the emergence of the non-integer occupancy 
eigenvalues shown in Figure [7| even when precisely two 
atoms are removed from the condensate state. 



VI. ENERGY OF TRANSLATIONAL MOTION 

In the previous sections we justified the separation of 
the quantum state into a center-of-mass translational mo- 
tion and relative motion of the atoms. We observed how 
the discrete excitations of the system in the case of strong 
attraction correspond to removal of one or more atoms 
from the condensate-like bound state of the system, while 
the lowest energy scale is associated with translation of 
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Figure 9: (Color online) Wave functions of the two particles 
excited out of the localized condensate state. For n = the 
wave function is real-valued and the two-particle wave func- 
tion amplitude is indicated by the color coding. 



the compound state. 

The localized state breaks the translational symmetry 
of the problem, while the restoration of the symmetry by 
a superposition of translated states gives rise to a simple 
analysis and interpretation of the ensuing energy. In con- 
tinuous space, Peierls and Yoccoz [25j, thus show that a 
superposition state with a complex phase variation with 
wave number K is a total momentum eigenstate and gives 
rise to the usual kinetic energy expression, fPK'^/2M, 
where M is the total mass of the system of particles. 
Peierls and Yoccoz also show that for an oriented quan- 
tum system, an angular superposition state weighted by 
Wigner rotation functions becomes an angular momen- 
tum eigenstate, e.g. of a nucleus or a molecule. In this 
case, however, the spectrum only recovers the rigid rotor 
form with terms proportional to L(L + 1) and M£ [1] 
if the state is sufficiently localized in Euler angles such 
that a Taylor expansion of the Wigner rotation functions 
is valid over the angular extent of the state. In nuclei, a 
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superfiuid core violates this, but the normal component 
may still be treated in the above sense and gives rise to 
rotational structure with an apparently reduced moment 
of inertia. 

In case of angular rotation of polarized electronic states 
in atoms, deviations from the rigid rotor energy structure 
in the form of an apparent increased moment of inertia, 
have been explained in [26j, by the angular extent of the 
states being too wide for the simple Taylor expansion of 
the Wigner function to be valid. 

The case of a finite, discrete lattice presents an interest- 
ing intermediate case between the continuous momentum 
and discrete angular momentum eigenstates, and the dis- 
persion relation for center-of-mass motion in the lattice 
is not a priori obvious. 

By our numerical diagonalization, we have already 
found the spectrum, but following [25j, we can also eval- 
uate the energy dispersion with direct reference to the 
representation of the translation and energy eigenstates 
I^J^) in terms of the underlying state given by (13). 
This leads to the expectation value of 



\H\^n) 



{<\Hfn 



lO 



(18) 



where we have used that commutes with H. Inserting 
the definition ( 14 ) into the numerator yields 



{a^\Hf„\a^) 



M-1 

M ^ 

k=0 



(a^lHT'^la^) (19) 



and a similarly expression with H replaced by 1 applies 
for the normalization factor in the denominator in ([l8|. 

Since the Hamiltonian favors nearest neighbor inter- 
actions, and since the states \a]^) become very narrow 
in the limit U/J — oo, we expect that only the terms 
k G {M — 1, 0, 1} contribute significantly in the sum (19). 
For the denominator we can make an even stronger argu- 
ment for the suppression of terms: Since \a^) is projec- 
tion onto the Fock basis states with center of mass in the 
wedge around the positive part of the real axis in the com- 
plex plane (see Figure [4|, and since T makes the center 
of mass rotate by the angle 2it/M in the complex plane, 
all terms {d^\T^\d^) with /c = 2, . . . , M — 2 vanish iden- 
tically, and terms with k = loik = M — l get their only 
contributions from Fock basis states components of 
with their center of mass phase angle exactly equal to 
±7r/M. We recall that this very strong truncation is valid 
because we chose to represent our states by Fock states 
with their center-of-mass coordinate within the wedge. 

Truncating (19)toA: = M — 1,0,1 and approximating 



the denominator in ( 18 ) by unity yields 



(20) 



which turns out to be in excellent agreement with the 
energy eigenvalues found from the numerical diagonal- 
ization of and our description thus accounts for the 
spectral features at all observable energy scales. 



VII. CONCLUSION 

In conclusion, we have solved the attractive Bose- 
Hubbard Hamiltonian numerically on a ring lattice with 
a finite number of atoms, and we have found that while 
the eigenstates have uniform mean population on all 
sites, the atomic wave function shows strong interpar- 
ticle correlations. These correlations are suggestive that 
the exact eigenstates can be represented as superposi- 
tion states of translated replica of compact (solitonic) 
composite many-body states. We have verified that this 
representation is, indeed, a useful one, and we have fur- 
ther demonstrated that the separation of center-of-mass 
and relative motion of the atoms explains features of the 
energy spectra very well. The argument is similar to 
the argument for the Born-Oppenheimer approximation 
applied in chemistry to separate nuclear and electronic 
motion. The separation is valid in the limit of strong 
interactions, where the energy scales become sufficiently 
different, and we observed how degeneracies in the system 
give rise to typical anti-crossing behavior with significant 
non-adiabatic coupling of different internal states. This 
phenomenon is also well known in Born-Oppenheimer 
theory, and plays a prominent role in the rich physics 
associated with conical intersections [27J. 

Our analysis supplements numerical calculations with 
a physical picture of the states of attractive bosons on a 
lattice. In this picture they form compact, macroscopi- 
cally populated condensate states, and they may be sub- 
ject to excitations of one or several atoms out of the con- 
densate state, which in turn occupies a superposition of 
different locations corresponding to center-of-mass quasi- 
momentum eigenstates of the system. 

Apart from offering insight in the anatomy of the 
many-body quantum state, the separation of the degrees 
of freedom in relative and center-of-mass motion paves 
the way for new mechanisms to excite and control the sys- 
tem dynamics. We imagine, for example, that the non- 
adiabatic coupling of different many-body states across 
a certain range of quasi- momentum eigenstates, may be 
applied to prepare interesting many-body superposition 
states. 



Appendix A: Ordering of Fock basis states 

In order to diagonalize a Hamiltonian which is given 
in a second quantization formulation as 



CstalcLt + ^ dstuvCLlalauCLy (Al) 

s,t s,t,u,v 



where the summation indices run over all the single par- 
ticle modes, we must write it as a matrix in a suitable 
basis, e.g., multi-mode Fock states, of many-body states: 
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Instead of iterating through ah pairs of basis states, a 
lot can be saved by using the sparsity of which only 
couples one basis state \Vt^) to a few others. 

For single particle operators we thus apply a\at to 
(for all s and t with non- vanishing Cst in ( |A1| )), and the 
resulting state can then be identified among the basis 
states. A similar procedure is used for the two-particle 
operators, but here one iterates over four indices instead 
of two. The identification procedure can be done conve- 
niently with a suitable hashing function that facilitates 
the handling of single mode indices in the full set of many- 
body states. Different general hashing functions are sug- 
gested in [28], but in the case where the Fock basis con- 
sists of all states with N bosons in M modes, [29j pro- 
poses a simpler hashing function that gives a convenient 
mapping between the single mode occupation numbers 
in the Fock basis states and the set of integers from 1 to 
. In [29] only the closed form formula for the hashing 
function is given, but here we define a binary relation to 
compare each two basis states and we derive the hash- 
ing formula from first principles. Finally, we show that, 
using the ordering induced by this particular binary re- 
lation, it is often not even necessary to use the hashing 
formula. 



1. Binary relation 

For two multi-mode A^-particle Fock states 

\Vt^) |mo,...,mM-i) 
\Vt^) = |no,...,nM-i) 

there will be a lowest k G {0,...,M — 1} such that 

ruj = Uj for j > k. 

We now introduce a binary relation ^at on the basis 

nk<mk \^a) 1 



rik > ruk 



(A2) 



This relation provides an ordering of the full set of basis 
states, and a simple recipe to find the successor of a given 
basis state, = |no, . . . , nM-i)- First find the highest 
A: G {0, . . . , M - 1} such that 

Uj = for j < /c, 

i.e. the state has the form 

= |0, . . . , 0, n/e, n/e+i, . . . , um-i) ' 

with n/g ^ 0. The next element in the ordered basis is 
then obtained by transferring one particle from the kth 
to the k + 1st mode, while the remaining particles in the 
kth mode are moved to the first one: 



This increment recipe clearly respects the ordering, 
l^a) l^a+i)' t)e proved inductively that 

successive use of this increment constructs all possible 
basis states. Thus, from now on we let 

B^ = {\n^),\n^),.-.,\^r^)) 

denote the ordered set of A/'-particle Fock basis states 



(A/' = 0, 1, 2, . . .) with the ordering induced by ( A2). 



2. Explicit enumeration of basis states 

In our calculations we make explicit use of the trans- 
lational invar iance of the problem, and hence we also 
need efficient means to identify the Fock state members 

of the equivalence classes {S^^'^}j_-^, and the action of 

the translation operator T. The action of T on a multi- 
mode Fock state is trivially carried out by shifting the 
occupation numbers cyclically, but the problem lies in 
identifying which basis state one arrives at. For this pur- 
pose it is advantageous to be able to calculate the basis 
state number directly (in the chosen ordering) without 
searching in a complete list of basis states. 

For this purpose, it is advantageous to represent 
the states by the A'-tuple (si, . . . , sat), where Sj G 
{0, . . . , M — 1} denotes which mode is occupied by the 
jth particle. This is a unique representation and it is 
not in conflict with the particle indistinguishability, if we 
choose si < S2 < - - - < sn- 

To each A^-tuple corresponds a basis state with a given 
number An{si^ . . . , sn) G {l, . . . , F])^} in our chosen or- 
dering, and to identify that number efficiently we note 
that the ordering begins with the states where we occupy 
only the first mode, then comes all possible ways of occu- 
pying only the first two modes, but with non-vanishing 
occupation of the second mode, then the first three modes 
and so forth. The last state which only occupies the first 
s modes 0, . . . , 5 — 1 is the state 



|0,...,0,A^,0,...,0) 



(A3) 



where the N appears in the place corresponding to mode 
5 — 1, and this state must have the number F^. The 
subsequent states, according to our ordering, start by 
putting one particle in the next mode and the remaining 
— 1 atoms in lower modes, and these states, all having 
one atom in the mode corresponding to the 5th mode, are 
to be numbered as FJ^ plus the sequential enumeration 
of the A" — 1 particles into the modes 0, . . . , s — 1: 

An{si, . . . , sn-i,s) = F^ + An-i{si, . . . , sn-i), 

where si < . . . < sn-i < s. By inserting this recursive 
relation into itself repeatedly we end up at the last term 



l^a+i) = - 1, . . . , 0, 0, n/c+i + 1, . . . , riM-i) • 



Ai(5i)=5i + l = Ff +1 
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(where we have extended the binomial coefficient such 
that (2) = when k > n > 0). This proves the expUcit 
number formula for the hashing function 



and since aj (restricted to Oj^) conserves the ordering of 



the Fock basis states in the sense (A5), we conclude that 
for ah |(^^^.)) e Of we have 



N 



(A4) 



Formula (A4) equips us with an explicit way of calcu- 



lating the number of the basis vector from its TV-tuple 
representation, and we can hence readily calculate the 
action of, e.g., the translation operator without having 
to search through the ordered list of states for a match 
of the shifted populations. 



Therefore, if we take two Fock basis states with occupa- 
tion of the jth and kth mode, respectively, G Of 

and \9.%^) G Of, we get 



3. Matrix elements 

The great advantage of using this ordering of the basis 
states is tha t, even though we just derived the explicit 
formula ( A4 ) , often we need not even calculate the hash- 



Ps -J Pt 

where the factor of proportionality is found from the oc- 
cupation numbers of the two states 



Ps Ps Pt Pt 



(A6) 



ing function. 

If we apply an annihilation operator dj to an N- 
particle Fock basis state |^^f ), with a non- vanishing oc- 
cupation in the jth mode, {ftf l^ijl^a ) 0? resulting 

state dj \nf) is proportional to an TV - 1-particle Fock ^he matrix representation of the operator d% is 



basis state, while dj \ ftf) vanishes if {^}f\nj\rtf) = 0. 
But the ordering has the property that if two A/'-particle 
states that are related as \^a) \^f) both occupy the 
jth mode, then the two N — 1-particle states obtained by 
operating with dj have the same ordering: 



-< 



N-l 



Nl 



(A5) 



thus found by identifying all the Fock basis states that 
have non-zero occupation of the jth and kth mode, re- 
spectively, while keeping their individual order. These 

( ■) 

states are indexed by the lists of numbers, {ps }s=i^ 
and {pi and then the only non- vanishing en- 

tries in the matrix representation of d^jdk are indexed as 



There is, indeed, a bijective correspondence between the 
ordered subset of basis states with non- vanishing occu- 
pation of the jth mode (j = 0, . . . , M — 1) 

of =f|nj,,),...,|nj,,),...,|nj,, }]cBff 
and the N — 1 particle multi-mode basis 



for s = 1, 

the numerical factors \AQ\ from the occupations of the 



, T](j_-\ , where we must calculate 



two modes in question for all the pairs of basis states. 



The ffist term in the Hamiltonian (Al) is merely a 



linear combination of such operators, so the matrix rep- 
resentation of that term follows easily from the above. 
The interaction terms with application of two annihila- 
tion operators and two creation operators benefit from 
the same use of ordered sets and convenient enumeration 
of states with non- vanishing occupancies. In the prepa- 
ration of the Hamiltonian matrix it is hence not even 
necessary to evaluate or utilize the full hashing function 
following from the preceding subsection. 
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